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ABSTRACT 


Turbulence is thought to be a primary driving force behind the early stages of 
star formation. In this framework large, self gravitating, turbulent clouds fragment into 
smaller clouds which in turn fragment into even smaller ones. At the end of this cascade 
we find the clouds which collapse into protostars. Following this process is extremely 
challenging numerically due to the large dynamical range, so in this paper we propose 
a semi analytic framework which is able to model star formation from the largest, giant 
molecular cloud (GMC) scale, to the final protostellar size scale. Due to the simplicity 
of the framework it is ideal for theoretical experimentation to explore the principal 
processes behind different aspects of star formation, at the cost of introducing strong 
assumptions about the collapse process. The basic version of the model discussed in 
this paper only contains turbulence, gravity and crude assumptions about feedback, 
nevertheless it can reproduce the observed core mass function (CMF) and provide 
the protostellar system mass function (PSMF), which shows a striking resemblance 
to the observed IMF, if a non-negligible fraction of gravitational energy goes into 
turbulence. Furthermore we find that to produce a universal IMF protostellar feedback 
must be taken into account otherwise the PSMF peak shows a strong dependence on 
the background temperature. 

Key words: stars: formation - turbulence - galaxies: evolution - galaxies: star for¬ 
mation - cosmology: theory 


1 INTRODUCTION 

Finding a comprehensive description of star formation has 
been one of the principal challenges of astrophysics for 
decades. Such a model would prove invaluable to under¬ 
standing the evolution of galactic structures, binary star 
systems and even the formation of planets. 

It has been long established that stars form from col¬ 
lapsed dense molecular clouds (McKee & Ostriker 2007a). 
Currently the most promising candidate for a driving process 
is turbulence, as it can create subregions with sufficiently 
high density so that they become self gravitating on their 
own, while also exhibiting close to scale free behavior (in ac¬ 
cordance with the observations of Larson 1981; Bolatto et al. 
2008). These fragments are inherently denser than their par¬ 
ents so they collapse faster, quasi independent from their 
surroundings. However, once they turn into stars they start 
heating up the surrounding gas (by radiation, solar winds 
or supernova explosions) preventing it from collapsing and 
forming stars (see Fig. 1). This process is inherently hierar¬ 
chical so it should be possible to derive a model that follows 
it from the scale of the largest self gravitating clouds, the 
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GMCs (~ 100pc), to the scale of protostars (~ 10“®pc). 
This is not possible in direct hydrodynamic simulations due 
to resolution limits, but can be treated approximately in 
analytic and semi-analytic models. 

This paradigm has been explored by Padoan et al. 
(1997) and Padoan & Nordlund (2002), then made more 
rigorous by Hennebelle & Chabrier (2008) who attempted 
to create an analytic model analogous to Press & Schechter 
(1974), which approximates the background density field as 
a Gaussian random field. A similar model was developed 
by Zamora-Aviles et al. (2012), however that did not rely 
on turbulence. Later Hopkins (2012a) expanded on these 
works by adopting the excursion set formalism to find the 
distribution of the largest self gravitating structures, which 
was found to be very similar to the observed distribution of 
GMCs. Similarily Hopkins (2012b) found that the distribu¬ 
tion of the smallest self gravitating structures fit well the 
observed CMF. Building on these results Hopkins (2013a) 
generalized the formalism to be applicable to systems with 
different equations of state and turbulent properties. 

Observed cores are sub-sonic and show no clear sign 
of fragmentation and the CMF looks very similar to 
the IMF apart from a factor of ~ 3 shift in the mass 
scale (Offner et al. 2014). However, if no other physics 
is assumed other than isothermal turbulence and grav- 
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Initial collapse of GMC 



Fragmentation of parent 


Independent collapse of fragments 



Figure 1. Evolution of collapsing clouds, with time increasing from left to right (darker subregions are higher-density, arrows denote 
regions which are independently self-gravitating and become thicker with increasing collapse rate). As the initial cloud collapses, density 
fluctuations increase (because gravitational energy pumps turbulence), creating self-gravitating subregions. These then collapse indepen¬ 
dently from the parent cloud, forming protostars at the end. These protostars can provide a sufficiently strong feedback that the rest of 
the cloud becomes unbound and ceases to collapse. 


ity, during the collapse the cores develop strong turbu¬ 
lence and eventually sub-fragment into smaller objects 
(Goodwin et al. 2004; Walch et al. 2012a, for discussion see 
Krumholz 2014). This implies that some additional physics 
must play a role, but there is no clear consensus on what 
it could be; magnetic fields (Nakano & Nakamura 1978; 
McKee & Ostriker 2007b), radiation (Krumholz 2011), cool¬ 
ing physics (Jappsen et al. 2005) etc. Using a cooling physics 
motivated “stiff” EOS Guszejnov & Hopkins (2015) incor¬ 
porated the time dependent collapse of the cores into the 
excursion set formalism and found that the distribution of 
protostars closely reproduced the observed IMF. 

These excursion set models did successfully reproduce 
the CMF, IMF and the GMC mass function, however they 
had several shortcomings. First, they did not account for the 
differences in formation and collapse times of clouds of differ¬ 
ent sizes (e.g. small clouds form faster and collapse faster). 
Secondly, the excursion set formalism describes the density 
field around a random Lagrangian point. This means that 
the spatial structure of a cloud can not be modeled directly 
(e.g. there is no way to find if a cloud forms binary stars). 
Finally, there is no self consistent excursion set model that 
follows from the GMC to the protostar scale (i.e. Hopkins 
2012 b covered scales between the galactic disk and cores, 
Guszejnov & Hopkins 2015 between cores and protostars). 
We believe these shortcomings can be overcome by moving 
away from the analytic excursion set formalism and instead 
adopting a simple semi-analytical approach with the same 
random field assumption. This framework would allow us 
to follow the evolution self gravitating clouds while resolv¬ 
ing both the GMC and protostellar scales and preserving 
spatial information. In this paper we will outline a possible 
candidate for such a model. 

The paper is organized as follows. Sec. 2 provides a gen¬ 
eral overview of the model, including the primary assump¬ 
tions and approximations and briefly outlines its numerical 
realization. Sec. 3 shows the simulated time evolution of the 


CMF and the protostellar system mass function (PSMF) 
which shows a striking similarity to the IMF. Sec. 3.2 also 
discusses the effects of having a temperature independent 
equation of state on the peak of the PSMF and the univer¬ 
sality of the IMF. Finally, Sec. 4 discusses the results and 
further applicability of the model. 


2 METHODOLOGY 

In short, instead of doing a detailed hydrodynamical sim¬ 
ulation involving gravity and radiation, our model assumes 
a simple stationary model for the density field, collapse of 
structures at constant virial parameter and an equation of 
state that depends on cloud properties. Starting from a 
GMC sized cloud it evolves the density field as the cloud 
collapses and pumps turbulence (this is not a bad approx- 
iation, see Robertson & Goldreich 2012; Murray & Chang 
2015; Murray et al. 2015). Note that our assumptions do not 
necessarily mean that all clouds have supersonic turbulence. 
Paper H has shown that if a medium has a “stiff” equation 
of state (7 > 4/3), then turbulence is dampened during col¬ 
lapse. Since it is observed that dense, low mass cores are 
subsonic while high mass, low density clouds are supersonic 
some form of physics is needed to remove the turbulent en¬ 
ergy. For that purpose we are using an equation of state 
that becomes stiff at high densities, which in combination 
with the constant virial parameter assumption makes dense 
clouds sub-sonic, arresting fragmentation. 

In the model, at each time step we search for self gravi¬ 
tating structures which we treat as new fragments, for which 
the process is repeated in recursion until a substructure is 
found that collapses to protostellar scale without fragment¬ 
ing. Our assumptions will be discussed in more detail in the 
following subsections while a step-by-step description of the 
algorithm is provided in Appendix A. 

Our model is a modified version of the excursion set 
model used by Guszejnov & Hopkins (2015) (henceforth re- 
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ferred to as Paper I) using the theoretical foundation of 
Hopkins (2013a) (henceforth referred to as Paper II). Due 
to the significant overlap between models we show only the 
essential equations and emphasize the differences and their 
consequences. If the reader is familiar with Paper I we sug¬ 
gest skipping to Sec. 2.3. 


2.1 The Density Field 

It is known that the density field in the cases of both sub 
and supersonic, isothermal flows follows approximately log¬ 
normal statistics (for corrections see Hopkins (2013b)). This 
means that if we introduce the density contrast (5(x) = 
In [p(x)/po] + S/2, with p(x) as the local density, po as the 
mean density and S' as the variance of Inp, it would follow 
a close to Gaussian distribution^, thus 

P(<5|S)«^exp(-y. (1) 

It is a property of normal and lognormal random variables 
that a linear functional of these variables will also be nor¬ 
mal/lognormal, thus the averaged density in a region has 
lognormal equilibrium statistics whose properties are pre¬ 
scribed by turbulence. Following Paper H this yields 

S(A)= [ AS(A)dlnA^ [ In [l-k (A)]dln A, (2) 
Jo Jo 

where A is the averaging scale, M (A) is the Mach number of 
the turbulent velocity dispersion on scale A and b is the frac¬ 
tion of the turbulent kinetic energy in compressive motions, 
which we take to be about 1/2 (this is appropriate for an 
equilibrium mix of driving modes, see Federrath et al. 2008 
for details. Paper I experimented with b ^ 1/4—1 and found 
no qualitative differences). 

It is important to note that although p is lognormal 
which means 5 is Gaussian, there is significant spatial corre¬ 
lation (i.e. p can not change instantly over arbitrarily small 
spatial intervals) so it is not possible to model the density 
field as a spatially independent random field. To circum¬ 
vent this issue we solve the problem in Fourier space since 
5 {k) is also lognormal, while there is little correlation be¬ 
tween modes so it is acceptable to assume them to be in¬ 
dependent (note: having correlated modes in Fourier space 
introduces only mild effects on the final mass functions, see 
Appendix A of Paper II for details). Combined with the 
fact that the number of modes in the [k^k dfc] range is 
dN(k) = (pi'Kk^Ak') Uk, where Uk is the mode density, we get 


^ It is a common misconception that analytical models such as 
the one presented in this paper take the total density distribution 
to be purely lognormal. While the density distribution in each 
cloud/fragment is indeed assumed to be locally lognormal on a 
single timestep, these have different means and deviations (see Eq. 
2) depending on their initial conditions and time, which means 
that the total distribution will be different. If we measure the 
density distribution in our calculations (see Fig. 2), we find it 
is approximately lognormal at low densities (set by the lowest 
density structure: the parent cloud), while the high mass end 
becomes a power law as it is a mass weighted average of the 
distributions for different substructures whose mass distribution 
is a power law (see Fig. 4). 



Figure 2. Time evolution of the distribution of density in a par¬ 
ent GMC of 10^ Mq. This is a mass weighted average of the den¬ 
sity distribution of all substructures in the parent cloud (which 
are all assumed to be lognormal with different parameters), thus 
the low mass end is set by the lowest density structure which is 
the parent cloud while the high mass end is a power law due to 
the power law like distribution of fragments (see Fig. 4). There 
is also a clear trend as the high mass end tail rises in time. This 
is caused by the formation of new self gravitating substructures 
(Federrath & Klessen 2013). 


the variance for an individual density contrast mode is 


5'mode(k) — 


ln{l + PMjkf) 


( 3 ) 


Paper H showed that to realize a steady state density 
contrast field with such variance and zero mean, the Fourier 
component 5(k, t) must evolve as 


5(k, t + At) = 5(k, t) (1 - At/Tk) + 

( 4 ) 

where 7?. is a Gaussian random number with zero mean and 
unit variance while Tk ~ Vt{k)/X is the turbulent crossing 
time on scale A ~ 1/fc, and the turbulence dispersion obeys 
Wi(A) oc A^“^ thus Tx oc A^^ (in our simulations we use 
p = 2, appropriate for supersonic turbulence, see (Murray 
1973; Schmidt et al. 2009; Federrath 2013)). 


2.1.1 The Equation of State 

It is easy to convince oneself that a purely isothermal or 
polytropic equation of state (EOS) would be a very poor de¬ 
scription of the complex physical processes contributing to 
the cooling and heating of clouds, however, modeling these 
processes in detail would require full numerical simulations. 
Instead we try to find a simple, heuristic EOS that captures 
the behaviors critical to our calculation. One of the most im¬ 
portant effect during collapse is the transitioning from the 
state where the cooling radiation efficiently escapes from the 
cloud to the state where the cloud becomes optically thick to 
it and heats up as it contracts. As the virial parameter is as¬ 
sumed to be constant, this leads to a decrease in turbulence, 
which effectively arrests fragmentation. This is essential to 
reproduce the IMF shape as pure isothermal collapse would 
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lead to an infinite fragmentation cascade. We adopt the same 
effective polytropic EOS model as Paper I where for small 
time steps (compared to the dynamical time): 


T{x,t + At) 


r(x,t) 


/ p{x,t +At) 
V p(x,t) 


7(t)-l 


( 5 ) 


where 7 {t) is the effective polytropic index of the cloud at 
time t. 

One of the main goals and advantages of our framework 
is that it allows the exploration of different physical EOS 
models simply and efficiently. For example, let us consider 
first the volume-density (n) dependent EOS model based on 
works like Masunaga & Inutsuka 2000; Glover & Mac Low 
2007 that follows the form 


1 

fO.8 

n < 

10 ® cm ® 


7(n) = 

1.0 

10 ® 

< < 10i° . 

(6) 

1 

[ 1.4 

n > 

10 ^° cm"® 



Simulations have shown that this leads to a ’turnover’ only 
at extremely low masses (~ 0.001 Mq, Fig. 9 later), making 
the IMF nearly a pure power-law at the observable masses. 
We will explore model and some of its physical consequences 
for observables in more detail in a future paper, but explic¬ 
itly show below that our semi-analytic model also captures 
this behavior. This is a valuable vindication both of the ac¬ 
curacy of the semi-analytic model (compared to full numer¬ 
ical simulations), and of the need for additional physics to 
establish the turn-over of the IMF. 

For purposes of this study, let us assume that we do 
not know the detailed origin of such physics (it may be due 
to magnetic fields, or radiative heating, for example, both 
of which we will explore in detail in follow up papers). The 
simplest approach, and one commonly adopted in numerical 
simulations, is to parametrize their effects via an ’effective 
equation of state’. Motivated by the work on radiative feed¬ 
back from [(Bate 2009; Krumholz 2011), let us consider a 
toy model where the effective EOS is not volume-density 
but surface-density (S) dependent: 


7(S) = 


0.7 


E < 3Mq/pc^ 

0.094In ( +0.7 3 < .. a < 5000 . 

\3MQ/pc-i y Mq/pc^ 

E > SOOOMq/pc^ 


1.4 


( 7 ) 


This is the same EOS as we used in Paper I. Note that the 
“turnover” where this becomes “stiff” is at much lower sur¬ 
face densities than we would obtain if we modeled cooling 
physics alone (Glover & Mac Low 2007) which would essen¬ 
tially give the same answer as our 7 (n) case above (for a com¬ 
parison of the two types of EOS models, see Guszejnov et al. 
(2015). Instead, we are assuming some form of physics makes 
the EOS stiffen at much higher surface densities - we choose 
the particular value here empirically, because it provides a 
reasonable fit to the observed IMF. We will then explore the 
consequences of such a parametrization, for the IMF and its 
time-evolution in different clouds 


2.2 Collapse: criterion and evolution 

It has been shown in Paper I and Paper II that the critical 
density for a (compared to the galactic disk) small, homoge¬ 
neous, spherical region of radius R to become self gravitating 


is 


PcTitjR) 

Po 


l+XeV 




+ .+^edge 



(8) 

where the two terms represent thermal and turbulent energy 
respectively. T{X) is the temperature averaged over the scale 
A, while To is the mean temperature of the whole collapsing 
cloud and we used the following scaling of the turbulent 
velocity dispersion and Mach number M 


M^{R) 


vKR) 

{cl (po)) 


= M; 


idge 



( 9 ) 


where Ro is the size of the self gravitating parent cloud and 
p is the turbulent spectra index, so the turbulent kinetic en¬ 
ergy scales as E{R) oc R^-, generally p € [5/3; 2], but in this 
paper, just like in Paper I we assume p = 2 as is appropriate 
for supersonic turbulence. 

It should be noted that the fragmentation process is 
complex even in the idealized case of homologous collapse 
(see Hanawa & Matsumoto 1999; Ntormousi & Hennebelle 
2015). This means that our method of finding self gravitating 
subregions using Eq. 8 is a strong approximation, however, 
a proper treatment would require drastically more computa¬ 
tion power which would go against one of the primary goal 
of the framework: the rapid exploration of parameter space 
and testing of physical models. 

Our goal is to create a model that resolves clouds from 
GMC to protostellar scales, so the initial structures of the 
model are the GMCs which themselves are self gravitating 
(first crossing scale in the excursion set formalism). This 
means they must satisfy Eq. 8 , which for spherical clouds 
{M{R) = ( 47 r/ 3 ) 77® p(i?)) in isothermal parents yields the 
mass-size relation: 


M = 




R 


2 


R 


sonic 



( 10 ) 


Note that for very high mass clouds a correction containing 
the angular frequency of the galactic disk would appear, 
however this term is small (see Paper II for details). Eq. 
10 introduces Tsonic which is the sonic length, the scale on 
which the turbulent velocity dispersion is equal to the sound 
speed, so in an isothermal cloud using the scaling of Eq. 9, 
we expect 

Rs.nic = ( 11 ) 

Meanwhile Maonic is defined as the minimum mass required 
for a sphere with 77sonic radius to start collapsing so 

2 r^R 

JWsomc — „ 

^coll ^ 

where G is the gravitational constant and Qcoii is the virial 
parameter for a sphere of the critical mass for collapse (see 
Eq. 15 later). For reasonable galactic parameters and tem¬ 
peratures i7sonic ~ 0.1 pc and Msonic ~ 6.5 Mq (assuming 
we use the value for Qcoii we specify in Sec. 2.2.1). 

Since the GMG in question has just started collapsing, 
the turbulent velocity at its edge must (initially) obey the 
turbulent power spectrum. Thus vl{R) oc R for the super¬ 
sonic and Vt{R) oc (the Kolmogorov scaling) for the 

subsonic case. Using the mass-size relation of Eq. 10 leads 
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to the following fitting function 

(1 + -Medge) -Medge ^ M 

1 + -^edge A^sonic 

which exhibits scalings of M oc Af® for the subsonic and 
M (X for the supersonic case respectively, and (cou¬ 
pled to the size-mass relation above) very closely repro¬ 
duces the observed linewidth-size relations (Larson 1981; 
Bolatto et al. 2008; Lada & Lada 2003). Note that dense re¬ 
gions will deviate from this scaling, as observed (see refer¬ 
ences above), because collapse “pumps” energy into turbu¬ 
lence (Robertson & Goldreich 2012; Murray & Chang 2015; 
Murray et al. 2015). 


2.2.1 Evolution of Collapsing Clouds 


One of the key assumptions of the previous models in Paper 

I and Paper II is that the kinetic energy of collapse pumps 
turbulence (Robertson & Goldreich 2012; Murray & Chang 
2015; Murray et al. 2015) whose energy is dissipated on a 
crossing time. As turbulent motion provides support against 
collapse, the collapse can only continue after this extra en¬ 
ergy has been dissipated by turbulence (see Sec. 9.2 in Paper 

II for details). This leads to the following equation for the 
contraction of the cloud: 


df 

df 



(14) 


where r{t) = R{t)/Ro is the relative size of the cloud at 
time t while r = t/to is time, normalized to the initial cloud 
dynamical time to ~ ‘^Q7oU^ {GMo/Rq) (see Paper II 
for derivation). In this case the initial dynamical time (to) 
and the crossing time only differ by a freely-defined order 
unity constant, so in our simulations we consider them to be 
equal without loss of generality. 

The other key assumption of the model is that collapse 
happens at constant virial parameter. We define Qcoii as 

Qcoll^ =4+ U? = (1 + M^dge) • (15) 

Note that Qcoii is not the Toomre Q parameter, merely the 
ratio of kinetic energy to potential energy needed to desta¬ 
bilize the cloud, thus the higher Qcoii the more unstable 
clouds are to fragmentation. One can find Qcoii using the 
Jeans criterion: 


0 > = [cl + Ut) ~ 4:7rGp, (16) 

which for the critical case (uj = 0) leads to 

One would be tempted to substitute in k — 2nfR, but 
that would be incorrect, as we have a spherical overden¬ 
sity with R radius to which the corresponding sinusoidal 
wavelength is not R. We therefore chose k ~ which 
yields Qcoii = 12/7r^ ~ 1.2. Note that all formulas contain 
Cs /Qcoii oc T /Qcoii so an uncertainty in the virial parameter 
is degenerate with an uncertainty in the initial temperature. 

Gombined, the above equations completely describe the 
collapse of a spherical cloud, as the EOS (Eq. 5-7) sets the 
temperature and thus the sound speed. Using that, Eq. 15 
provides the edge Mach number, which allows us using Eq. 
14 to calculate the contraction speed. 


2.3 Differences from previous models 

So far we are following the same assumptions as Paper I 
and Paper II, however, instead of simulating a stochastic 
density field averaged on different scales around a random 
Lagrangian point (the basis of analytic excursion set models) 
we use a grid in space and time. This means that we directly 
evolve the 5 (k) modes to simulate the density field. This 
allows us to preserve spatial information as we now have 
information about the relative positions and velocities of 
substructures. 

Having a proper density field not only allows us to take 
basic geometrical effects into account (as substructures are 
still assumed to be spherical) but it allows a proper applica¬ 
tion of the self gravitation condition of Eq. 8. The excursion 
set formalism finds the smallest self gravitating structure a 
point is embedded in. The problem is that this “last crossing” 
structure may have further self gravitating fragments which 
do not contain the aforementioned point. These substruc¬ 
tures will form protostars of their own (see Fig. 1) leaving 
their parent cloud with less mass which in turn might not be 
self gravitating anymore. This is not addressed in excursion 
set models which instead simply assume 100% of the mass 
endiug up in protostars of different sizes (which of course is 
not realistic), while the proposed grid model predicts only 
about 5% (see Sec. 3. 2), which in fact depends on the physi¬ 
cal assumptions of the model (i.e. how to deal with unbound 
material). 

It should be noted that like the model of Paper I, in 
this first study we include no explicit feedback mechanism. 
Instead the model utilizes a few crude approximations to 
account for the qualitative effects of feedback. First, it is 
assumed that the clouds that becomes unbound by frag¬ 
mentation stop collapsing and “linger” for a few dynami¬ 
cal times (during which they may form new self gravitating 
fragments) before being heated up/blown up/disrupted by 
feedback from the newly created protostars in such a fashion 
that they can no longer participate in star formation^. Note 
that this assumption is made for convenience, it is not in¬ 
herent in the code as it is possible to implement direct feed¬ 
back prescriptions. Similarly magnetic fields are neglected 
in this base model, but can be easily implemented into the 
framework. Like in Paper I we neglected the effects of ac¬ 
cretion and protostellar fragmentation when comparing to 
the IMF as the protostellar system mass function (from now 
on PSMF) is already a good enough qualitative fit so their 
effects are assumed to be modest (except for the very high 
and low mass ends where fragmentation could provide a high 
mass cut off while accretion could affect the turnover point, 
see McKee & Offner 2010 for details on the protostellar mass 
function). We would also like to note that it is possible to 
apply a crude implementation of supernova feedback by sim¬ 
ply stopping the evolution after a few Myrs (when enough 
supernovae have exploded to unbind the GMG). Since the 
simulation provides a time dependent output, it can be done 
during post-processing. Of course, the point of our frame- 


^ For example photoionization can destroy the molecular 
cloud (Dale et al. 2012; Walch et al. 2012b; Geen et al. 2015), 
while both supernovae (Iffrig & Hennebelle 2015) and outflows 
(Arce et al. 2007) can provide momentum for turbulence or eject 
material. 
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work is that one could easily add models for feedback, and/or 
accretion if desired. 

We would like to note that using hydrodynamical sim¬ 
ulations would allow a much more realistic treatment of 
certain details of the problem, however the large dynamic 
range (10“® — 100 pc) and the long range gravitational in¬ 
teractions make such attempts extremely computationally 
intensive, preventing one from getting substantial statistics. 
A further issue with direct hydrodynamical simulations is 
that they involve the full, detailed form of all physical in¬ 
teractions, making it harder to pinpoint the primary driving 
mechanisms behind certain phenomena. 

In summary we propose a semi-analytical model which 
has negligible computational cost but still captures phenom¬ 
ena (e.g. spatial correlation, motion of objects, complicated 
time dependence) which are beyond the capabilities of the 
analytical excursion set formalism. Our intention in this pa¬ 
per is not to present a “complete” model of star formation, 
but rather to illustrate the power of this approach with a 
first study involving only turbulence and self-gravity. 


3 EVOLUTION OF THE IMF AND CMF IN 
GMCS 

In this section we present an application of the model 
for simulating the collapse of an ensemble of GMCs (dis¬ 
tributed following the first crossing mass function obtained 
by Hopkins 2012b, see Fig. 3). This includes simulating a 
number of GMCs of different masses where the initial con¬ 
ditions are set by Eq. 9 and Eq. 10. The clouds are assumed 
to start with fully formed turbulence (as GMCs form out of 
an already turbulent medium) which means that before sim¬ 
ulating the collapse the density field is initialized to have the 
appropriate lognormal distribution. The output of the code 
contains the formation time and properties (e.g. mass, posi¬ 
tion, velocity) of individual protostars along with snapshots 
of the hierarchical structure of bound objects at different 
times. In Sec. 3.1 we investigate the latter and compare the 
distribution of nonfragmented structures with the observed 
CMF. Later, in Sec. 3.2 we discuss the time evolution of 
PSMF and how it relates to the IMF and whether it can be 
universal without invoking feedback physics. 

3.1 Fragmentation and self-gravitating 
substructures: the observed CMF 

It is well known that during their collapse clouds fragment 
into smaller self gravitating structures (see Fig. 1). It is in¬ 
structive to see how much mass is bound in structures of dif¬ 
ferent sizes. Fig. 4 shows the time evolution of the number of 
structures of different sizes counting all “clouds-in-clouds”, 
which follows a distribution similar to the observed IMF and 
CMF (for quick overview see Offner et al. 2014), however it 
has a significantly shallower slope® of roughly The 

distribution is established fairly quickly and is maintained 

® In this paper the approximate high mass end behavior is esti¬ 
mated by fitting a power law between O.SMq-IOOMq. The error 
presented in the figures only account for the uncertainty in the 
fitting. 



Figure 3. Initial mass function of GMCs according to the excur¬ 
sion set model of Hopkins (2012b) compared to the observations 
(X symbols) and empirical fitting function (dashed black line) of 
Rosolowsky (2005). The normalization of the plot is arbitrary. 



Figure 4. Time evolution of number of bound structures of dif¬ 
ferent masses in a parent GMC of 10® Mq. Here we count all 
self-gravitating structures, including clouds embedded in other 
clouds, cores etc. The plot is normalized so that integrated mass 
(f M j corresponds to the mass of gas bound in self 

gravitating clouds relative to the total mass of the parent GMC, 
which explains the decreasing trend with time as more and more 
gas ends up in either protostars or becomes unbound. The upper 
end cuts off close to the parent GMC mass. The high mass power 
law fitting is done according to Footnote 3. 


until the collapse of the parent cloud ends. This mass func¬ 
tion of bound structures is consistent with the cloud in cloud 
picture shown in Fig. 1 in that there is a vast hierarchy of 
bound structures embedded in each other. 

Observationally finding the substructure of a GMC is 
very challenging (although see Rosolowsky et al. 2008), most 
observers instead concentrate on the so called cores which 
are collapsing clouds that have no self gravitating fragments. 
Figure 5 shows the total CMF (time and mass averaged 
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Figure 5. Comparison of the average simulated CMF with the 
observed CMF by Sadavoy et al. (2010) in different clouds in the 
Milky Way (the plot is normalized so that the peak of the CMF 
is set to unity). Note that observations which are below the com¬ 
pleteness limit are also included (see the original paper for de¬ 
tails). The simulated CMFs are averaged both over time (assum¬ 
ing the age of GMCs is uniformly distributed in the [0,5] Myr 
range) and the CMC mass function (following Fig. 3). The dif¬ 
ferent initial critical masses in this case reflect having different 
T/Qcoll values, for definition see Eq. 20. 


over an ensemble of GMCs following the distribution shown 
in Fig. 3) for different inital parameters. The simulated 
CMF reproduces the shape of observed results, having both 
a turnover point and a slightly shallower high mass slope 
(~ than the canonical Salpeter result of ~ 

for the IMF (see Offner et al. 2014). 

Fig. 6 clearly shows that there is very small differ¬ 
ence between the CMF turnover masses and high mass 
slopes between GMCs of different sizes after 1 Myr. This 
is because early collapse is roughly isothermal so these 
clouds all have the same characteristic fragment mass (Merit, 
see Eq. 20 for details). Systems which are on the same 
linewidth-size relation (i.e. they form out of the same tur¬ 
bulent cascade) will always have the same Maonic, Merit (see 
Hennebelle & Chabrier 2008; Hopkins 2012b). During later 
evolution the GMCs heat up at a different pace as the dy¬ 
namical times are different. Meanwhile Fig. 7 shows that 
there is a clear trend of increasing turnover mass with time 
in each cloud. This phenomenon and its possible cause is fur¬ 
ther investigated in Sec. 3.2. This trend is not visible in case 
of the physical EOS of Eq. 6 as the peak is well below the 
stellar mass scales (see Eig. 9). Nevertheless, this scenario 
shows that in the absence of a dominant Merit the initial 
CMF turns over around the sonic mass scale (as shown by 
previous analytical works e.g. Hennebelle & Chabrier 2008; 
Hopkins 2012b), but this mass scale gets “forgotten” during 
the fragmentation cascade. 


3.2 Evolution of the PSMF 

We now examine the mass function of the final collapsed 
objects, the protostellar system mass function (PSMF). 



Figure 6. The CMF in GMCs of different masses 1 Myr after col¬ 
lapse starts for each cloud (using EOS of Eq. 7). The plot is nor¬ 
malized so that integrated mass corresponds to the relative mass 
of gas bound in cores, the peaks are denoted with solid circles. 
The high mass power law fitting is done according to Footnote 
3. Both the turnover mass and the high mass slope exhibit very 
little sensitivity to the mass of the parent GMC similar to what 
was found by Hennebelle & Chabrier (2008); Hennebelle (2012). 


In Fig. 8 we show that parent clouds of all masses pro¬ 
duce similar to Salpeter scalings the high mass end with 
lower mass clouds producing slightly steeper slopes. Also, 
there is a clear trend of increasing turnover mass with in¬ 
creasing parent mass, unlike the case of the CMF (See Fig. 
6 ). It is worth noting that the GMC mass function is top 
heavy, which means that the high mass clouds dominate 
the integrated mass function. If we accept this result then 
it suggests a possible observational bias of the IMF as most 
observations focus on smaller clouds in the Milky Way. Also, 
turbulent fragmentation does not produce a cloud mass de¬ 
pendent “maximum stellar mass”. 

The increasing turnover mass for both PSMF and CMF 
is related to the equation of state. In a turbulent cloud, self 
gravitating fragments of different sizes form, which (accord¬ 
ing to the EOS of Eq. 7) have different effective polytropic 
indices. According to the EOS there exists a threshold in the 
surface density (Ecrit) above which 7 > 4/3, stabilizing the 
cloud against further fragmentation. Thus it is instructive to 
find the critical mass (Merit) corresponding to Ecrit. Using 
the collapse condition of Eq. 8 and expanding up to linear 
order in 7 around 1 (this is a good approximation during 
most of the cloud’s lifetime as the collapse starts at close to 
isothermal conditions) yields that E > Ecrit requires that 

7 

R < i?crit = Ro - 7 -r-, (18) 

^(l + -^eV)--^eV+7-l 

where R is the fragment radius and i?o, Eq, 7 = 7 (Eq) are 
the radius, surface density and the effective polytropic index 
of the parent cloud. From Eq. 18 we can find the critical 
mass Merit = Ecrit below which fragments are unlikely 
to collapse (note; according to the EOS of Eq. 7 the critical 
surface density Ecrit ~ 2400 Mq/pc^). These formulas can 
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Figure 7. Left: Time evolution of the CMF in a 10® Mq parent GMC using the 7 (S) EOS of Eq. 7. The plot is normalized so that 
integrated mass corresponds to the mass of gas bound in self gravitating clouds relative to the total mass of the parent GMC, which 
explains the downwards trend since less and less gas is bound in cores as more protostars are produced and the cloud gets heated by 
contraction. The high mass power law fitting is done according to Footnote 3. There is a clear trend in the turnover mass (the peaks are 
denoted with solid circles) which increases significantly while preserving the overall shape of the function (e.g. high mass slope). Right: 
Time evolution of the CMF in a 10“^ Mq parent GMC using the physically motivated EOS of Eq. 6 (a density dependent EOS where the 
transition point to the 7 > 1 regime is calculated from cooling physics). As expected the CMF has a peak around the sonic mass at early 
times, however, that feature gets “washed out” by the fragmentation cascade which is not arrested by this EOS until very small scales. 



Figure 8 . Protostellar system mass function (PSMF) after col¬ 
lapse ends in parents of different masses assuming our simple 
equation of state. The Salpeter slope is always present (the high 
mass power law fitting is done according to Footnote 3). For these 
assumptions there appears to be “too many” brown dwarfs, and 
too much dependence on the parent GMC mass. These are the 
direct consequences of the EOS of the gas. 


be simplified by assuming isothermal collapse (7 1 ) and 

that the parent GMC is highly supersonic (Algdge ^ 1)5 Eq. 
11 yields then: 


r> -^oEq So 

-ttcrit ~ 1 ^ — -tX-sonic ^ 

^edge^crit E, 


'crit 


(19) 


Using the mass-size relation of Eq. 10 and that Rq ';$> iZsonic 


we obtain 

Ecrit 167ri?2„„i,Ecrit 47rG2 Ecrit’ 

( 20 ) 

The critical mass only depends on the cloud temperature 
and the equation of state. A similar sensitivity to the initial 
temperature has been found by Bate (2009) using a Jeans 
mass argument. Assuming that there exists a critical density 
Pcrit where some physics terminates the fragmentation cas¬ 
cade the corresponding Jeans mass will simply be oc 
It is easy to see that this is the same result one would get 
when trying to find the critical mass using a 'y{n) EOS. 

Fig. 9 shows the time evolution of the time and ensemble 
averaged PSMF for different initial Merit values (the differ¬ 
ent critical masses in these cases arise from having different 
cr/QcoiiEcrit; where we fix Qcoii and Ecrit and vary Tinit, 
for definition see Eq. 20) which all produce a shape similar 
to the IMF but with different peak masses. If we compare 
the results to the canonical IMF fitting functions of Kroupa 
(2002) and Chabrier (2005), then it is clear that the average 
PSMF always reproduces the Salpeter scalings however the 
turnover point is heavily influenced by T/QcoiiEcrit . Since 
Qcoii is a constant this implies that the average temperature 
of the cloud could have a signiheant effect on the turnover 
point if Ecrit is constant. Meanwhile, Fig.9 also shows that 
the physical EOS of Eq. 6 has such a low characteristic mass 
that the resulting PSMF in the stellar mass range is just a 
power law. Nevertheless, the position of the peak is still sen¬ 
sitive to the initial conditions (oc if one extends the 

plot to substellar mass scales. 

Fig. 10 shows how this critical mass evolves in time for 
our default model assumptions (Ecrit = const.). It is clear 
that Merit correlates well with the peaks of the PSMF of the 
corresponding time interval. 

This increase of the critical mass with time has an in¬ 
teresting consequence. Fig. 11 shows that the average time 
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Figure 9. Evolution of the averaged PSMF (normalized to in¬ 
tegrated mass) for different initial critical masses (set by having 
different T'/Qcoil^crit values, for definition see Eq. 20) compared 
to results using the “traditional” EOS of Eq. 6 and the canonical 
IMF of Kroupa (2002) and Chabrier (2005). The PSMF is aver¬ 
aged both over time (assuming the age of GMCs is uniformly dis¬ 
tributed in the [0,5] Myr range) and the GMC mass function (fol¬ 
lowing Fig. 3). We included the standard Merit = 0-03 Mq (solid 
red), an Merit = 0.08 Mq (solid blue) and an Merit = 0.2 Mq 
(solid black) scenarios with the 7 (S) EOS along with a run which 
had the physically motivated 7 (n) EOS of Eq. 6 . For realistic tem¬ 
peratures (10 — 30 K) the critical mass of the latter is well below 
the stellar mass range so the PSMF becomes a pure power law. 
Meanwhile, for the 7 (E) EOS case the PSMF shape is similar for 
different critical masses, and there is a clear shift of the peak to 
higher masses with increasing Merit- Iri all cases the high mass 
end is close to the Salpeter result. 



Figure 10. The peak masses of the PSMF of different time inter¬ 
vals (solid line with symbols) and the critical mass (dashed lines) 
for different parent GMC masses according to Eq. 18. The criti¬ 
cal mass correctly predicts the qualitative evolution of the peak 
mass. 


Figure 11. Average time of formation for protostars of differ¬ 
ent masses (the error bars represent the standard deviation) in 
a model with an invariant EOS. There is a clear trend of more 
massive protostars forming at later times (which is consistent with 
the shifting of the turnover mass in Fig. 10), however the scatter 
is comparable to this difference. Nevertheless it is clear that most 
massive stars only start forming after roughly a Myr after the 
cloud starts collapsing. Changing this requires additional physics 
beyond turbulence, gravity and cooling. 


of formation monotonically increases with the protostellar 
system mass. 

So, if the equation of state does not depend on tem¬ 
perature (e.g. our 7 (E) is invariant) then the turnover mass 
shows a strong (oc T^) dependence on the initial conditions 
which would likely lead to a non-universal IMF (oc in 
the 7 (n) case). A possible solution to this issue is if Ecrit from 
Eq. 20 has a temperature dependence. This perfectly plausi¬ 
ble, just recall that the effective EOS is just a crude approx¬ 
imation of complex cooling physics. Bate (2009) argues that 
radiative feedback effectively weakens the dependence of the 
Jeans mass on density, making the turnover mass less sensi¬ 
tive to initial conditions. A similar example is provided by 
Krumholz (2011), where the initially formed protostar ’seed’ 
heats up its environment, preventing it from collapsing. This 
dense cloud is heated up to Theatmg oc ^ by 

the accretion luminosity from the protostar^, which, using 
our EOS language, roughly translates to Ecrit oc which 
would produce a constant Merit, and thus a universal IMF. 

In a paper in preparation we will explore this feedback 
model in a fully spatially-dependent framework. For now, let 
us consider a simple experiment where Ecrit cx T^. 

Fig. 12 compares the results of two simulations, one with 
Ecrit = const, and one with Ecrit cx T^. Although the latter 
still shows some time dependence, the shifting of the peak 
is greatly reduced, making it more consistent with observa¬ 
tions, even though the only assumption about feedback was 
that it prevents collapsed cores from accreting from their 


^ One can derive this temperature by assuming an optically thick 
cloud in equilibrium that is heated by accretion luminosity Lace ~ 
MT ~ M/tffT oc and cooled by thermal radiation 

icooi ~ 
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surroundings. Note that our aim with this experiment was 
only to demonstrate what would be required from a purely 
EOS based model to produce an invariant IMF, any other 
physics that sets the critical mass of the EOS constant would 
achieve similar results. 

An important question of star formation is what frac¬ 
tion of the gas ends up in stars. The analytical excursion 
set models like in Paper I could not answer that question 
as they assume by default that 100 % of the mass ends 
up in bound structures similar to the Press-Schecter model 
(Press & Schechter 1974) of dark matter halos which they 
are based on). However, our semi-analytic framework here 
allows us to explore different assumptions for the time- 
dependent behavior of both bound and unbound gas, and 
thus (in principle) to make predictions for this quantity. 

In the “basic” models presented in this paper, we as¬ 
sume that whenever a core collapses and forms a star, any 
remaining mass in its parent cloud which is no longer self- 
gravitating (once the core is fully collapsed) is simply thrown 
out of the system. This is meant to represent a very crude toy 
model for the effects of feedback (from e.g. protostellar jets) 
on the parent sub-clumps from which the stars form. With 
this assumption, we find an integrated star-formation effi¬ 
ciency (after all mass either turns into stars or is unbound) 
of ~ 5 — 10% for GMCs of all sizes. Interestingly, this is al¬ 
most completely independent of the EOS we assume (either 
constant Ecrit or Ecrit oc H), as long as it terminates the 
fragmentation cascade at roughly the same point. Of course, 
if we assume this gas remains bound to the total system, so 
it is simply recycled back to the “top level” of the original 
fragmentation hierarchy until it is consumed (which obvi¬ 
ously corresponds to a no-feedback case), then we trivially 
predict that eventually all gas turns into stars. Of course, 
the effects of realistic feedback are much more complex than 
these simplistic assumptions, and we could adopt arbitrar¬ 
ily complex models (for example, evolving each protostar 
and tracking explicitly location-dependent photo-ionization 
feedback, which we then use to explicitly calculate whether 
gas is unbound from the system). We note this result simply 
to demonstrate the utility of these semi-analytic models for 
rapidly exploring different assumptions regarding the effects 
of feedback. 


4 CONCLUSIONS 

The aim of this paper is to provide a general framework 
for the modeling of star formation through turbulent frag¬ 
mentation from the scale of GMCs to the scale of stars in 
order to quickly test the effects of different assumptions 
and new physics. Such a tool could allow theorists to ex¬ 
plore different models and parameters before committing 
significant resources towards a detailed numerical simula¬ 
tion. We propose a semi analytical extension of the model of 
Guszejnov & Hopkins (2015) (Paper I) that we believe is de¬ 
tailed enough to capture the physics essential for modeling 
the formation of stars without being too demanding numer¬ 
ically. Just like the analytical excursion set models it does 
not simulate turbulence directly, instead it assumes that the 
density follows a locally random held distribution whose pa¬ 
rameters evolve in time so that virial equilibrium is satished. 
This is an assumption about turbulent collapse that needs 


to be tested in future work. The density held is directly re¬ 
solved on a grid which preserves spatial and time information 
allowing the implementation of more detailed physics (e.g. 
proper checking for self gravitation, time dependent cloud 
collapse) and the analysis of the spatial structure. This is 
not possible in the excursion set formalism which describes 
the density held around a random Lagrangian point. This 
also means that unlike the analytical models not 100 % of 
the mass ends up in protostars. 

The presented form of the model contains only the min¬ 
imally required physics (turbulence, self gravity, some equa¬ 
tion of state). It is however possible to integrate more so¬ 
phisticated models to provide a more accurate description of 
these processes. Also, since the output of our model contain 
the time dependent evolution of the CMF and the PSMF, 
one can easily apply corrections during post processing to ac¬ 
count for ehects like protostellar fragmentation or supernova 
feedback (stop the evolution when enough SNe exploded). 

By applying this framework to modeling the collapse of 
giant molecular clouds, we found that even the basic model 
qualitatively reproduces the observed core mass function. 
The GMF evolution has little dependence on the mass of 
the parent GMC mass. 

Another result of the simulation is the mass distribution 
of all bound structures in the cloud. This appears to have the 
same shape as the GMF with a shallower slope of roughly 
M~^'^ at the massive end. These clearly show the hierarchy 
of bound structures. 

One of the main results of our basic model is the pro¬ 
tostellar system mass function (PSMF) which is obtained 
by following the collapse of an ensemble of GMCs following 
a GMC mass function determined by Hopkins (2012b). As 
in Paper I we found that the PSMF is qualitatively very 
similar to the observed IMF: it exhibits a close to Salpeter 
slope almost independent of the initial conditions, while the 
turnover mass is mainly set by the equation of state and the 
initial temperature. 

Due to the minimalistic nature of the model we man¬ 
aged to pinpoint the physical quantities influencing the dif¬ 
ferent features of the PSMF and thus the IMF. We found 
that the Salpeter slope at the high mass end is a clear con¬ 
sequence of turbulence (as shown before in Paper I) where 
the inclusion of extra physics only causes slight deviation 
from the pure power law behavior. Furthermore we found 
that in a medium with a stiff equation of state the actual 
turnover point in leading order is set by the local tempera¬ 
ture (Merit OC T^/Ecrit). 

We found that if we assume a 7 (E) equation of state 
then the PSMF for protostars of the same age changes as 
the parent cloud collapses: the turnover mass increases with 
time. This can be explained by the increase of Merit. This 
leads to a quadratic dependence of the turnover mass on the 
initial temperature which is inconsistent with the observed 
universality of the IMF. This means that it is not possible 
to derive a universal IMF with an equation of state that 
has no temperature dependence. One way to ’fix’ the model 
is by implementing the feedback from protostars. Using the 
assumptions of Krumholz (2011) in leading order the heat¬ 
ing from the protostars cancel the aforementioned quadratic 
scaling (due to Ecrit oc T^), leading to a close to universal 
turnover mass. 
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Figure 12. PSMF for protostars in a parent GMC of 10® Mq for an EOS with S^rit = const, (left) and for an EOS with Sj-rit oc 
(right). The solid circles show the peaks, which move considerably less for the Ecrit oc case. As implied by Eq. 20, if S^rit oc then 
^crit ~ const, and the IME becomes invariant. 
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APPENDIX A: BASIC SIMULATION 
ALGORITHM 

In this appendix we detail step-by-step how the basic version 
of the simulation works (see flowchart of Fig. Al), but note 
that it can be greatly expanded with new physics, as long as 
the fundamental assumption (locally random density modes) 
is kept. 

(1) We begin with a GMC sized cloud whose initial parame¬ 
ters (mass, radius, temperature, density, edge Mach number, 
sound speed etc.) are derived from its mass (M), the sonic 
mass (Msonic) and length (Rsonic), using the mass-size rela¬ 
tion of Eq. 10 and linewidth-size relation of Eq. 9. These are 
all initialized on a 3D spatial grid, of resolution N x N x N 
chosen such that the final statistics converge (we found this 
happens at A > 16). The density field is initialized assum¬ 
ing that it is lognormal (variance set according to Eq. 3) 
using the full density power spectrum model (transforming 
to Fourier space and back), while the temperature field fol¬ 
lows the density according to the desired equation of state. 

(2) We take timestep At {At tdyn and At <C tcross(d), where 
d = 2R/N is the spatial resolution of the grid). This means: 

(a) Global contraction of the cloud (all scales shrink, density 
uniformly increases) according to Eq. 14. 

(b) The density perturbation power spectrum (5(k) is up¬ 
dated following Eq. 4, which assumes density mode statis¬ 
tics obey a local “random walk” in phase space. The actual 
density field is calculated by Fourier transforming to real 
space and normalizing the field with the cloud mass (this 
way mass is conserved). 

(c) The temperature field is updated according to new den¬ 
sities and the chosen EOS (see Eq. 5). 

(d) The cloud scale Mach number is updated according to 
our assumption that the virial parameter is constant dur¬ 
ing collapse. 

(3) We now check whether any self-gravitating substructures 
have formed by using a Monte Carlo method that involves 
placing spheres of all possible sizes at random positions and 
testing them using the collapse criterion of Eq. 8. 

(4) If such a region is found it is “removed” temporarily and 
expanded into its own grid. This new grid will have a higher 
spatial resolution than its parent, thus density modes on the 
newly available small scales need to be initialized (larger 
modes are inherited from previous grid). We then repeat 
steps (2)-(4) on this new grid. This means that during the 
evolution of its fragments the parent cloud is “frozen” in 
time. This is motivated by the fact that the dynamical time 
of fragments is smaller as tdyn oc 1/ y/p, so they evolve “fast” 
compared to their parents. Note that all clouds keep track of 
physical time, so it is possible to properly date the formation 
times of protostars and clouds. 

(5) The time evolution of each cloud/grid continues until: 

(a) The cloud reaches the protostellar size scale {R < Rmin), 
below which it is assumed to have formed a protostar. 

(b) The cloud is still self gravitating after a number of dy¬ 
namical times {t > tmax)®. After this limit is reached the 
cloud is assumed to have cooled and collapsed through 

® This can happen if 7 > 1, as f in Eq. 14 does not reach zero in 
a finite amount of time. 


other means. Essentially, this represents non-fragmenting 
cores. 

(c) The cloud stops being self-gravitating. This can happen 
if a cloud loses enough of its mass that it becomes un¬ 
bound. Since virial equilibrium is enforced this means no 
turbulence, which means no more fragmentation. In the 
model presented above these clouds are not forming stars 
or contributing to the mass of the protostars forming from 
their fragments, instead this material is “thrown away” 
(this represents “feedback” in some sense, see Sec. 2.3). 
Note that it is possible within the framework to return 
this unbound material to the parent GMC where it may 
form stars, but for simplicity in the presented model we 
chose not to do that. 

(6) Clouds that formed protostars are removed the properties 
of the protostars are cataloged. We then return to the parent 
cloud and continue its evolution from Step (4). 

(7) This continues until 100% of the original original mass of 
the cloud is either in protostars or unbound. The final output 
is the catalog of protostars. Note that it is also possible to 
get the CMF by exporting the properties of bound structures 
at a specified time. The whole process is repeated for large 
number of initial GMCs (with different random seeds) to 
gain adequate statistics. 
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Figure Al. Basic algorithm of fragmentation code. The bold numbers in each box show which step from Appendix A they represent. 
See Appendix A for more detailed description. 
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